
#on crick
R
#list.dirs<-function(path){
#	x<-dir(path) 
#	x[file_test("-d",paste(path,x,sep=""))]
#}

##function to get call rate###
callrate <- function(x){ 
	#round(sum(x<0.05,na.rm=T)/sum(!is.na(x)),6)
	round(sum(x<0.05,na.rm=T)/length(x),6)
}


#path<-"/nt/ressrv08programs3/BMI/zxs/research/methylation/tcga/BRCA/450k/lvl3/"
#outpath<-"/nt/ressrv08programs3/BMI/zxs/research/methylation/tcga/BRCA/450k/lvl3/merge"

path<-"/data2/bsi/RandD/Arraybased_RND/METHYLATION/tcga_zs/BRCA/450k/lvl3/"
outpath<-"/data2/bsi/RandD/Arraybased_RND/METHYLATION/tcga_zs/BRCA/450k/lvl3/merge"


#dirs<-c("jhu-usc.edu_BRCA.HumanMethylation27.Level_1.1.0.0","jhu-usc.edu_BRCA.HumanMethylation27.Level_1.2.0.0","jhu-usc.edu_BRCA.HumanMethylation27.Level_1.3.0.0","jhu-usc.edu_BRCA.HumanMethylation27.Level_1.4.0.0")
dirs<-dir(path)
dirs<-dirs[grep("Level_3",dirs)]

##################merge file by columns in directory ##################
for (d in 1:length(dirs)){
	samplefiles<-list.files(paste(path,dirs[d],sep=""))
	samplefiles<-samplefiles[grep("BRCA",samplefiles)] ## skip any other files
	cat("directory:", dirs[d],"\n")
	for (i in 1:length(samplefiles)){
	#for (i in 1:2){ ##test for two files
		cat("directory:", dirs[d], "and file:", samplefiles[i],"\n") 
		if (d==1 & i==1){
		methy<-read.table(paste(path,dirs[d],"/",samplefiles[i],sep=""),header = T, skip = 1, sep = "\t",stringsAsFactors = F)
		#methy$beta<-round(methy[,2]/(methy[,2]+methy[,5]),4)
		#methy<-methy[,c(1,13,2,5,12)]
		methy<-methy[,1:2] #detection p value does not have much meaning all 0
		colnames(methy)<-paste(strsplit(samplefiles[i],split="\\.")[[1]][6],c("CpG","AvgBeta"),sep=".")
		#escape \\ for .
		}else{
		tmp <- read.table(paste(path,dirs[d],"/",samplefiles[i],sep=""),header = T, skip = 1, sep = "\t",stringsAsFactors = F)
		#tmp$beta<-round(tmp[,2]/(tmp[,2]+tmp[,5]),4)
		tmp<-tmp[,1:2]
		colnames(tmp)<-paste(strsplit(samplefiles[i],split="\\.")[[1]][6],c("CpG","AvgBeta"),sep=".")
		methy<-merge(methy,tmp,by.x=colnames(methy)[1],by.y=colnames(tmp)[1])
		}
	}
}
colnames(methy)
dim(methy) #  485577    596  -- 347 samples; only has average beta
methy[1:4,1:3]
names(methy)[1]<-"CpGID" 
names(methy)<-sub(".AvgBeta","",names(methy))
#average beta
#avgbeta<-methy[,c(1,seq(2,dim(methy)[2],4))]
#names(avgbeta)<-sub(".Avg_Beta","",names(avgbeta))
#pvals<-methy[,c(1,seq(5,dim(methy)[2],4))]
#png(paste(outpath,"avgbeta.boxplot.png",sep=""),width=2400,height=960)
#boxplot(avgbeta[,-1],xlab="samples",ylab="beta",cex=0.1)
#dev.off()

#wite out for later data analysis
#write.table(methy,paste(outpath,"methydata.fourcolms.each.txt",sep=""),sep="\t",quote=F, row.names=F)
RowNAs<-apply(methy[,-1],1,FUN = function(x){sum(is.na(x))})
table(RowNAs)
#  0    595
#449042  36535 #36535 CpGs are marked as NA in all samples
ColNAs<-apply(methy[,-1],2,FUN = function(x){sum(is.na(x))})
table(ColNAs)
methy<-na.omit(methy)
dim(methy) # 449042    596
png(paste(outpath,"/avgbeta.boxplot.png",sep=""),width=2400,height=960)
boxplot(methy[,-1],xlab="samples",ylab="beta",cex=0.1)
dev.off()

write.table(methy,paste(outpath,"/merged.avgbeta.txt",sep=""),row.names=F,sep="\t",quote=F)
methy<-read.table(paste(outpath,"/merged.avgbeta.txt",sep=""),sep="\t",header=T,stringsAsFactors=F,check.names=F,row.names=1)
#barcode<-strsplit(names(methy[-1]),"-")
#samtype<-unlist(lapply(barcode, FUN=function(x) return(x[4])))
#samtype<-sub("A","",samtype)
#samtype<-sub("B","",samtype)
#unique(samtype) #"01" "11" "06" 01-09: tumor; 10-19: normal; 20-29: control samples
#samtype.str<-ifelse(samtype %in% c("01","06"),"tumor","normal")
#methytest<-t(methy)
saminfo<-read.table(paste(outpath,"/BRCA.mappings.use.txt",sep=""),sep="\t",header=T,stringsAsFactors=F,check.names=F)

###get tumor/normal pairs ####
#tnpairs<-unlist(subset(saminfo,WithNormalPair==TRUE,TCGAID))
tnpairs<-subset(saminfo,!is.na(Pair))
dim(tnpairs) #181   8
#sample TCGA-E2-A15K has two tumors and one normal which makes 181 so we have 90 tumor/normal pairs
#TCGA-E2-A15K-01A-11D-A12R-05	TCGA-E2-A15K
#TCGA-E2-A15K-06A-11D-A12R-05	TCGA-E2-A15K
#TCGA-E2-A15K-11A-13D-A12R-05	TCGA-E2-A15K
tnpair.beta<-methy[,names(methy) %in% tnpairs$TCGA.Barcode]
write.table(tnpair.beta,paste(outpath,"/merged.avgbeta.tnpairs.txt",sep=""),sep="\t",quote=F)

###can start from here to save time###
#test<-read.table(paste(outpath,"/merged.avgbeta.tnpairs.txt",sep=""),sep="\t",header=T,check.names=F)
tnpair.beta<-read.table(paste(outpath,"/merged.avgbeta.tnpairs.txt",sep=""),sep="\t",header=T,check.names=F)
##transpose the data
#methy<-as.matrix(methy)
#methy<-t(methy)
#methy.saminfo<-merge(sameinfo,methy,by.x="TCGA.Barcode",by.y="row.names")

##use CpGassoc package#####
library(CpGassoc)
tnpair.beta[1:4,]
tnpairs[1:4,]
table(tnpairs$TCGA.BATCH)
table(tnpairs$histology)
tnpairs$TisType<-ifelse(tnpairs$histology=="Normal Tissue","Normal","Cancer")

#differentially methylatied CpGs between tumor normal, no covariate adjusterment
tndiffCpGs<-cpg.assoc(tnpair.beta,factor(tnpairs$TisType,levels=c("Normal","Cancer"),ordered=T),large.data = T)
#summary(tndiffCpGs)
#plot(tndiffCpGs)
str(tndiffCpGs) #5 lists
tndiffCpGs$results[1:5,]

####annotations####
#intalled in my windows but would not install in crick###
#source("http://bioconductor.org/biocLite.R")
#biocLite("IlluminaHumanMethylation450k.db")
#ls("package:IlluminaHumanMethylation450k.db")
#x <- IlluminaHumanMethylation450kCHR37
#xx <- as.list(x)[1:10]
#xx
#read in annotation
####focus on chrX only#####
anno450k<-read.table(paste(outpath,"/GPL13534_HumanMethylation450_15017482_v.1.1.KeyColms.txt",sep=""),sep="\t",header=T,check.names=F,stringsAsFactors=F,quote = "")
anno450k.chrX.IlmnID<-subset(anno450k,CHR=="X")[,1]
length(anno450k.chrX.IlmnID) #11232
tnpair.beta.chrX<-tnpair.beta[which(row.names(tnpair.beta) %in% anno450k.chrX.IlmnID),] ###some CpG were filtered out by TCGA
dim(tnpair.beta.chrX) #10852   181
write.table(tnpair.beta.chrX,paste(outpath,"/merged.avgbeta.tnpairs.chrX.txt",sep=""),sep="\t",quote=F)






